function [ correctedImageSet alpha beta scoreMatrix negValues magVariableCorrection]...
    = computeRedImageBgCorrection(imageSet, imageFeatures)
%computeRedImageBgCorrection
%   Computes a two variable optimized correction of two retinal images based on the red
%   reflection or IR image. Optimization is based on minimizing the variance
%   of the ratio image. NOTE: negValues occur when scoreMatrix = 100
%
%
%   Inputs
%   ------------
%
%   imageSet must be a vector that contains the images used in the algorithm.
%   It contains three images defined as imageSet.blueIm,
%   imageSet.otherIm, and imageSet.redIm. imageSet.blueIm will be the
%   numerator in the ratio, imageSet.otherIm will be the denominator in the
%   ratio, and imageSet.redIm will be the image used for background
%   correction.
%
%   imageFeatures is a vector that contains information needed in our algorithm: 
%   the fovea center, the macular region of interest, and the contrast regions for
%   background correction. The imageFeatures vector contains individual
%   vectors defined as imageFeatures.fovRow, imageFeatures.fovCol,
%   imageFeatures.iMacRegion, and imageFeatures.iContrastRegion.
%
%   Outputs
%   ------------
%
%   correctedImageSet is a vector that contains the output of corrected images
%   This vector has three images defined as correctedImageSet.blueIm,
%   correctedImageSet.otherIm, and correctedImageSet.redIm.
%
%   alpha is the percentage of the red image used to correct the background
%   in the blue image. alpha is chosen by minimizing the variance of the
%   ratio image, which is the blueCorrected/otherCorrected image.
%
%   beta is the percentage of the red image used to correct the background
%   in the yellow image. beta is chosen by minimizing the variance of the
%   ratio image, which is the blueCorrected/otherCorrected image.
%
%   bestScore is the variance of the image that optimized alpha and beta.
%
%   scoreMatrix contains all the variance calculations of the ratio
%   images generated from the 2-D grid search of alpha and beta.
%
%   Description
%   ------------- 
%   computeRedImageBgCorrection is used to isolate the image signal of
%   macular pigment in the eye so that the pixel intensities correspond
%   to the relative concentrations of macular pigment. PLEASE ADD MORE
%   DESCRIPTION WITH MORE ACCURACY.
%
%   Example
%   ------------
%   Read in blue, other, and red image
%   imageSet.blueIm = imread('blue.tif')
%   imageSet.otherIm = imread('other.tif')
%   imageSet.redIm = imread('red.tif')
%
%   Define the contrast regions used for correction by user input. See
%   MpMapperGUI code for more detail on how these regions are defined.
%   imageFeatures.iContrastRegion = iContrastRegion
%
%   Compute 2-D variable background correction
%   [correctedImageSet alpha beta bestScore scoreMatrix] = ...
%     computeRedImageBgCorrection(imageSet, imageFeatures)


% Below is the X range optimization
% handles.A_RANGE  = [0.01 0.6]; % MULTIPLE OF REF. IMAGE SUBTRACTED FROM BLUE AV.
% handles.B_RANGE  = [0.01 0.6]; % MULTIPLE OF REF. IMAGE SUBTRACTED FROM YEL  AV.
% handles.DELTA    = 0.01;       % INCREMENT 
% 
% handles.A_vals  = handles.A_RANGE(1) : handles.DELTA : handles.A_RANGE(2);
% handles.B_vals  = handles.B_RANGE(1) : handles.DELTA : handles.B_RANGE(2);
% 
% handles.X_vals  = zeros(length(handles.A_vals)*length(handles.B_vals),1);
% count = 1;
% for i = 1 : length(handles.A_vals)
%     for j = 1 : length(handles.B_vals)
%         A = handles.A_vals(i); B = handles.B_vals(j);
% 
%         X = (A - B)/(1 - B);        
% 
%         handles.X_vals(count) = X;
% 
%         count = count + 1;
%     end
% end
% 
% handles.X_RANGE = [min(handles.X_vals) max(handles.X_vals)];
% 
% handles.X_vals = handles.X_RANGE(1) : handles.DELTA : handles.X_RANGE(2);
% 
% best_X = 0;
% best_score = 10^6;
% 
% %This for loop will determine which x value is good by saying x is just
% %below 10^6
% for X = handles.X_vals
%     
%     image_difference = handles.blueAvgIm(handles.iMacBoundary) ...
%         - X*handles.redIm(handles.iMacBoundary) - (1-X)*handles.otherAvgIm(handles.iMacBoundary);
%     
%     score = mean(abs(image_difference(:)));
%         
%     if (score < best_score)
%         best_score = score;
%         best_X = X;
%     end
%         
% end
% 
% disp(' ');
% disp('BEST X VALUE:')
% disp(['X = ' num2str(best_X) ', where X = A - B*[(1-A)/(1-B)].']);
% disp(['score  = ' num2str(best_score) ...
%         ', where score = MEAN( N(blue) - X*N(red) - (1-X)*N(yellow) ) over macular boundary.']);
% 
% % Given optimal value for X (above), select optimal B (and associated A)
% % to minimize variation in [N(blue) - A*N(red)]/[N(yel) - B*N(red)] 
% % over selected contrast region.
% 
% best_A = 0;
% best_B = 0;
% best_ratio_score = 10^6;
% score = zeros(length(handles.A_vals), length(handles.B_vals));
% i = 1;
% j = 1;
% for a = handles.A_vals
%     for b = handles.B_vals
%         handles.blue_corrected = handles.blueNormalizedIm(handles.iContrastRegion) - a.*handles.redIm_normalized(handles.iContrastRegion);
%         handles.yel_corrected = handles.yellowNormalizedIm(handles.iContrastRegion) - b.*handles.redIm_normalized(handles.iContrastRegion);
%         
%         ratio_im = handles.blue_corrected ./ handles.yel_corrected;
%         if any(ratio_im < 0)
%             ratio_score = 100;
%         else
%             ratio_score = var(ratio_im);
%         end
%         
%         score(i,j) = ratio_score;
%         j = j+1;
%         if ratio_score < best_ratio_score
%             best_ratio_score = ratio_score;
%             best_A = a;
%             best_B = b;
%         end
%     end
%     i = i+1;
%     j = 1;
% end



alphaRange = [0.01 0.6]; % Range of alpha values in 2-D grid search
betaRange = [0.01 0.6]; % Range of beta values in 2-D grid search
delta = 0.01; %Step of values in 2-D grid search

alphaVals = alphaRange(1) : delta : alphaRange(2);
betaVals = betaRange(1) : delta : betaRange(2);

bestAlpha = 0;
bestBeta = 0;
bestRatioScore = 10^6; % initialize to worst possible variance value
scoreMatrix = zeros(length(alphaVals), length(betaVals));
i = 1;
j = 1;
rowRange = imageFeatures.fovRow - imageFeatures.macRadius: imageFeatures.fovRow + imageFeatures.macRadius;
colRange = imageFeatures.fovCol - imageFeatures.macRadius: imageFeatures.fovCol + imageFeatures.macRadius;
for a = alphaVals
    for b = betaVals
% Below code is for optimizing along contrast region. This sets the region
% outside of the 4 degrees equal to the average value of the pixels along
% the 4 degree perimeter.
        blueCorrected = imageSet.blueIm - a.*imageSet.redIm;
        blueCorrected(imageFeatures.iNonMacRegion) = mean(blueCorrected(imageFeatures.iMacBoundary));
        otherCorrected = imageSet.otherIm - b.*imageSet.redIm;
        otherCorrected(imageFeatures.iNonMacRegion) = mean(otherCorrected(imageFeatures.iMacBoundary));
% Compute the ratio image and set the region outside 4 degrees equal to the
% pixel values along the 4 degree perimeter in the ratio image. This is
% done to so that we do not account for variations outside the 4 degree
% region.
        ratioIm = blueCorrected ./ otherCorrected;
        ratioIm(imageFeatures.iNonMacRegion) = ...
            mean(ratioIm(imageFeatures.iMacBoundary));
% Determine in their are negative pixel values. If so we need to remove
% these results since they are implausible. We will remove by setting their
% total variation to a high number (rememeber we are optimizing bto the
% lowest total variation).
        negValueLocations = any(ratioIm < 0);
        numNegValues = sum(negValueLocations);
        if numNegValues > 10
            ratioScore = 100000;
        else
            ratioScore = compute_total_variation(ratioIm(rowRange,colRange)); 
        end
        
        scoreMatrix(i,j) = ratioScore;
        j = j+1;
        if ratioScore < bestRatioScore
            bestRatioScore = ratioScore;
            bestAlpha = a;
            bestBeta = b;
        end
    end
    i = i+1;
    j = 1;
end

alpha = bestAlpha;
beta = bestBeta;
bestScore = bestRatioScore;

correctedImageSet.blueIm = imageSet.blueIm  - alpha*imageSet.redIm;
correctedImageSet.otherIm  = imageSet.otherIm  - beta*imageSet.redIm;

% Renormalize the corrected images according to Dr. Bonner's method
correctedImageSet.blueIm = correctedImageSet.blueIm.*(1/(1-alpha));
correctedImageSet.otherIm = correctedImageSet.otherIm.*(1/(1-beta));

% Calculate the number of negative scores in the 2-D grid serach
negValues = sum(sum(scoreMatrix == 100));


%Generate a box and compute the average value at the center of the fovea
boxsize = 10;
rowRange = imageFeatures.fovRow - boxsize:imageFeatures.fovRow + boxsize;
colRange = imageFeatures.fovCol - boxsize:imageFeatures.fovCol + boxsize;
pixelValues = double(imageSet.redIm(rowRange, colRange));

%Compute magnitude of correction at center and boudnary
magVariableCorrection.blueCenter = mean2(alpha.*pixelValues);
magVariableCorrection.otherCenter = mean2(beta.*pixelValues);
magVariableCorrection.blueBoundary = mean(alpha.*...
    imageSet.redIm(imageFeatures.iMacBoundary));
magVariableCorrection.otherBoundary = mean(beta.*...
    imageSet.redIm(imageFeatures.iMacBoundary));








